#include <stdio.h>
#include "spline_lib.h"
#include "randomlib.h"
#include "step_down.h"

/* global data */
int iknot;
struct L2_1D_DATA *  dat;
struct SP_1D * sp, * save_sp[1000];
int  num_rand_splines;

int main(int argc, char** argv)
{
  FILE *fp, *fp1;

  int i, k, order, n, nobs, inform, n_large, n_final, off_set, order_vec[500], rank[500];
  int seed1, seed2, count, sum, i_maic, i_true_ssq, num_step_dn, num_fn;

  struct SP_1D  * sp1;
  struct L2_1D_DATA *  dat_actual;

  double h, x, y, sse, cv, cor, glob_sse=999999, glob_sse_loc=999999., sigma, act_err, aic, bic;
  double mod_aic, mod_aic_vec[500], act_err_vec[500], sse_vec[500], dim_vec[500], sort_v[500];
  double min_maic, min_true_ssq;
  float read_sigma;

  fp = fopen("step_output.txt", "w");
  fp1 = fopen("step_input.txt", "r");

  fscanf (fp1, " %d %d %f",  &seed1, &seed2, &read_sigma);
  fscanf (fp1, " %d %d %d",  &num_fn, &nobs, &order);
  fscanf (fp1, "%d %d %d", &count, &num_step_dn, &n_large);
  
  /* internal parameters for the approximant is 
     order, num_step_dn, count, num_step_dn, n_large 
     which should be part of the input parameters
  */

  sigma = read_sigma;

  h = 1./(nobs -1);
 
  /*  initialize the data struct  */
  dat = data_1d_initialize(  nobs, &inform) ;
  dat_actual = data_1d_initialize(  nobs, &inform) ;

  /* initialize the random numbers */
  RandomInitialise(seed1, seed2);

  /* generate random 2-d data */
  for(i=0; i<dat->nobs; i++){
    dat->weight[i] = 1;
    dat->xdata[i]  = i*h;
    x = dat->xdata[i];
    dat->ydata[i]  = fun(x, num_fn) + sigma*RandomGaussian(0., 1.);

    dat_actual->weight[i] = 1;
    dat_actual->xdata[i] = i*h;
    dat_actual->ydata[i] = fun(x, num_fn);
  }

  k= -1;

  /* The main loop to do the approximation */
  /* parallelizable */
  for(order = 2; order <7; order++)
    {

      for(n_final = order; n_final < count+order; n_final++)
	{

	  k++;
	  sp1 = step_down_iter(num_step_dn,  order, n_large, n_final);
	  sse = get_L2_error( sp1, dat);
	  
	  save_sp[k] = spline_1d_copy(sp1);
	  act_err = get_L2_error( sp1, dat_actual);

	  aic = dat->nobs*log(pow(sse,2)/dat->nobs) + 2*(2*sp1->n - sp1->order);
	  bic = dat->nobs*log(pow(sse,2)/dat->nobs) + log(dat->nobs)*(2*sp1->n - sp1->order);

	  mod_aic = aic + (2*sp1->n - sp1->order);
	  mod_aic_vec[k] = mod_aic;
	  act_err_vec[k] = act_err*act_err/dat->nobs;

	  sse_vec[k] = sse*sse/dat->nobs;
	  order_vec[k] = sp1->order;
	  free_SP_1D(sp1);
	}
    }


  i_maic = 0;
  i_true_ssq = 0;
  min_maic = mod_aic_vec[0];
  min_true_ssq = act_err_vec[0];

  for(i=1; i< count*5; i++){
    if(min_maic > mod_aic_vec[i]){
      i_maic = i;
      min_maic =  mod_aic_vec[i];}
    if(min_true_ssq > act_err_vec[i]){
      i_true_ssq = i;
      min_true_ssq =  act_err_vec[i];}
  }


  for(i=0; i< count*5; i++)
    sort_v[i] = act_err_vec[i];

  gnomesort(5*count, sort_v);

  k = 0;
  for(i=0; i<5*count; i++){
    if (sort_v[i] == act_err_vec[i_maic]) break;
    k = k++;}


  /* write to the log file if required */

  fprintf(fp, "M_aic spline has rank = %d \n \n", k);
  fprintf(fp, "\n \n The rn parameters are: \n");
  fprintf(fp, "seed1 = %d, seed2 = %d, sigma =  %f \n \n", seed1, seed2, sigma);
  fprintf(fp, "Function number = %d  and the value at 0 is  %f \n", num_fn, fun(0., num_fn));
  fprintf(fp, "NOBS = %d \n \n", dat->nobs);
  fprintf(fp, "order = %d \n \n", sp1->order);
  fprintf(fp, "num_step_dn = %d ", num_step_dn);
  fprintf(fp, "n_large = %d \n \n \n", n_large);

  if(i_true_ssq == i_maic) 
    fprintf(fp, "M_aic correcly predicts min IMS \n \n");
  else
    fprintf(fp, "No agreement between best fit and aic \n \n");

  fprintf(fp, "act_err/nobs  mod_aic      SSE/nobs  \n");

  i = 0;

  for(order = 2; order <7; order++){
    for(k= 0; k< count; k++){
      fprintf(fp, "%e     %f   %f  %d\n", act_err_vec[i], mod_aic_vec[i], sse_vec[i], order_vec[i]);
      i++;}	
    fprintf(fp, "\n");
  }


sum = 0;
i = 0;
for(order = 2; order <7; order++){
	for(k=0; k<count-2; k++){
	if(sse_vec[i+1] > sse_vec[i]) 
	{sum++; break;}
	i++;}
	if (sum > 0 ) break;
	i = count*(order-1);

}
	
if(sum >0) {
		fprintf(fp, "                           Caution errors NOT ordered  \n");
		fprintf(fp, "                           order = %d   position = %d  %d \n", order, k, i);
			}
	else 
		fprintf(fp, "                           Errors ordered as expected \n");

fprintf(fp, "\n \n");
fprintf(fp, "act_err/nobs  mod_aic      SSE/nobs    Spline order\n \n");
k = sp->order;
i = i_true_ssq;
fprintf(fp, "Actual Best Error \n");
fprintf(fp, "%f     %f      %f         %d \n \n", act_err_vec[i], mod_aic_vec[i], sse_vec[i], order_vec[i]);
i = i_maic;
fprintf(fp, "Minimum Modified AIC \n");
fprintf(fp, "%f     %f      %f         %d \n \n \n", act_err_vec[i], mod_aic_vec[i], sse_vec[i], order_vec[i]);


fprintf(fp, "\n \nThe M_AIC spline  values are:\n");
for(i=0; i<dat->nobs; i++){
	x = dat->xdata[i];
	y = sp_1d_value(save_sp[i_maic], dat->xdata[i], 0);
	fprintf(fp, "%f %f %f\n", x, y , dat->ydata[i]);
}

/*
fprintf(fp, "\n \n");
for(k=0; k< count; k++){
fprintf(fp, "\n the %d -th splines is \n", k);
for(i=0; i<dat->nobs; i++){
	x = dat->xdata[i];
	y = sp_1d_value(save_sp[k], dat->xdata[i], 0);
	fprintf(fp, "%f %f %f\n", x, y , dat->ydata[i]);
}
}
*/

 return 0;
} 


void get_random_knots(struct SP_1D * sp, struct L2_1D_DATA * data ){
int i, num_int_knots;
double *  temp, low, high;

num_int_knots = sp->n - sp->order;
temp = (sp->knot + sp->order);

/*  Note this assumes that the xdata field is ordered */
low = data->xdata[0];
high = data->xdata[data->nobs -1];

for(i=0; i<num_int_knots; i++)
	temp[i] = low +RandomUniform()*(high - low);
	
gnomesort(num_int_knots, temp);

/*
for(i=0; i<num_int_knots; i++)
	sp->knot[sp->order + i] = temp[i];

free(temp);
*/
}


void gnomesort(int n, double * ar){
	int i = 0;
	while(i<n){
		if (i==0 || ar[i-1] <= ar[i]) i++;
		else {double tmp = ar[i]; ar[i] = ar[i-1]; ar[--i] = tmp;}
	}

}

double get_L2_error(struct SP_1D * sp, struct L2_1D_DATA * data ){
int i;
double sse;	
sse = 0.;
for(i=0; i<data->nobs; i++)
sse += (data->weight[i])*pow(data->ydata[i] - sp_1d_value(sp, data->xdata[i], 0), 2);

return sqrt(sse);
}

double opt_knot(void){
	int   j;
	double sse_old, sse, a, b, tol; 
	
	tol = .0001;
	
sse_old = get_L2_error( sp, dat );
for(j=0; j<99; j++){
	for(iknot = sp->order; iknot < sp->n  ; iknot++){
		a = sp->knot[iknot-1], b = sp->knot[iknot+1];
		sp->knot[iknot] = fminbr(a, b, f, tol);
	
	}

	for(iknot = sp->n -1; iknot >= sp->order  ; iknot--){
		a = sp->knot[iknot-1], b = sp->knot[iknot+1];
		sp->knot[iknot] = fminbr(a, b, f, tol);
		
	}
	

	for(iknot = sp->n -1; iknot >= sp->order  ; iknot--){
		a = sp->knot[iknot-1], b = sp->knot[iknot+1];
		sp->knot[iknot] = fminbr(a, b, f, tol);
		
	}

		for(iknot = sp->order; iknot < sp->n  ; iknot++){
		a = sp->knot[iknot-1], b = sp->knot[iknot+1];
		sp->knot[iknot] = fminbr(a, b, f, tol);
	
	}
		sse = get_L2_error( sp, dat );
	/*	printf( " sse = %f   sse_old = %f \n", sse, sse_old);  */
		if( fabs( (sse-sse_old)/(sse_old + .0001) ) < .0001 )
			break;
		sse_old = sse;
}

	sse = get_L2_error( sp, dat );
	return sse;
}

double f(double t){
	double sse;
	sp->knot[iknot] = t;
	get_L2_1D_spline(dat, sp);
	sse = get_L2_error( sp, dat );
	return sse;
}

	

struct SP_1D * par_fit(int seed, int order, int n, double * glob_sse){
int i,j, k, inform;

struct SP_1D  * sp1, *sp2;

double   sse, glob_sse_loc=999999.;

/* Data and sp are global structures */

*glob_sse=999999;


k = seed;

/*
printf(" \n \n The seeds for this job are %d and 15000\n", k);
printf("The number of random splines is %d per iteration \n", num_rand_splines);
printf("Order = %d and dimension = %d \n\n ", order, n);

*/

sp = sp_1d_initialize( order, n,  &inform);

sp_1d_set_knots(sp, 0., 1.);

get_L2_1D_spline( dat, sp );
sp1 = spline_1d_copy(sp);
sp2 = spline_1d_copy(sp);


RandomInitialise(k,15000);

/* Begin major loop  */

for(k=0; k<9; k++){

	glob_sse_loc=999999.;

for(j=0; j< num_rand_splines; j++){
get_random_knots(sp, dat);



get_L2_1D_spline( dat, sp );

sse = get_L2_error( sp, dat);





if(sse < glob_sse_loc){
	free_SP_1D(sp2);
	glob_sse_loc = sse;
	sp2 = spline_1d_copy(sp);}
		
/*	  printf("The error is %f \n", sse);  */

}


/*
printf("The iteration is %d \n", k);

printf("The global best error before opt is %f \n", glob_sse_loc);
for(i=0; i< sp2->n + sp2->order; i++)
	printf("%f ", sp2->knot[i]);

printf("\n ");

*/

free_SP_1D(sp);
sp = spline_1d_copy(sp2);

sse = opt_knot(); 

if(sse < *glob_sse){
	free_SP_1D(sp1);
	*glob_sse = sse;
	sp1 = spline_1d_copy(sp);
	
/*
	printf("The global best error after opt is %f \n", *glob_sse);
for(i=0; i< sp1->n + sp1->order; i++)
	printf("%f ", sp1->knot[i]);


*/


	}

printf("\n \n");

}

return   sp1;
}




struct SP_1D * step_down_fit( int order, int n_large, int n_final, double * glob_sse){
int i,j, k, isave, inform;

struct SP_1D  * sp1, *sp2;

double   sse, glob_sse_loc=999999.;

/* Data and sp are global structures */

/*
  printf(" \n \n The seeds for this job are %d and 15000\n", k);
  printf("The number of random splines is %d per iteration \n", num_rand_splines);
  printf("Order = %d and dimension = %d \n\n ", order, n);
*/

sp = sp_1d_initialize( order, n_large,  &inform);

sp_1d_set_knots(sp, 0., 1.);

get_L2_1D_spline( dat, sp );
sp1 = spline_1d_copy(sp);
sp2 = spline_1d_copy(sp);
get_random_knots(sp,  dat );


/* Begin major loop  */

k = n_large;

while( k > n_final){

	glob_sse_loc=999999.;
	
/*	if(5*(k/5) == k)
		sse = opt_knot(); 
*/
for(i = order; i<k; i++){

	free_SP_1D(sp2);
	sp2 = get_one_less_knot(i, sp);	
    get_L2_1D_spline( dat, sp2 );

    sse = get_L2_error( sp2, dat);


if(sse < glob_sse_loc){
	isave = i;
	free_SP_1D(sp1);
	glob_sse_loc = sse;
	sp1 = spline_1d_copy(sp2);}
/*	  printf("The error is %f \n", sse);  */

}
 //printf("For k = %d  The error is %f \n", k, glob_sse_loc);
 k--;
 free_SP_1D(sp);
 sp = spline_1d_copy(sp1);

} //end while
/*
printf("The iteration is %d \n", k);

printf("The global best error before opt is %f \n", glob_sse_loc);
for(i=0; i< sp2->n + sp2->order; i++)
	printf("%f ", sp2->knot[i]);

printf("\n ");

*/

free_SP_1D(sp);
sp = spline_1d_copy(sp1);

sse = opt_knot(); 

if(sse < *glob_sse){
	free_SP_1D(sp1);
	*glob_sse = sse;
	sp1 = spline_1d_copy(sp);	}
	
/*
	printf("The global best error after opt is %f \n", *glob_sse);
for(i=0; i< sp1->n + sp1->order; i++)
	printf("%f ", sp1->knot[i]);


*/




//printf("\n \n");
return   sp1;

}


struct SP_1D * get_one_less_knot(int i, struct SP_1D * sp){
struct SP_1D * sp1;
int j, inform;
	
	sp1 = sp_1d_initialize( sp->order, sp->n -1,  &inform);

	for(j=0; j< sp->n + sp->order; j++){
		if(j < i)
			sp1->knot[j] = sp->knot[j];
		
        else if(j>i)
			sp1->knot[j-1] = sp->knot[j];
		
	}

	return sp1;	
}




struct SP_1D * step_down_iter(int num_tries, int order, int n_large, int n_final){
int i, j;
double glob_sse = 999999, cur_sse;
struct SP_1D * sp_final, * sp1;

cur_sse = glob_sse;
for(i=0; i< num_tries; i++){
sp1 = step_down_fit(order, n_large, n_final, &glob_sse);
if(glob_sse < cur_sse){
	sp_final = spline_1d_copy(sp1); 
	cur_sse = glob_sse;} 

}
free_SP_1D(sp1);
return sp_final;
}


double fun(double x, int i){
	double y= 0.;
	switch(i){
		case 1: y = 10*(sin(x*x*x*12  ));  break;
		case 2: y = 1/(.1 + sin(x*x*8 )*sin(x*x*8 )); break;
		case 3: y = fabs(cos(x*x*6)); break;
		case 4: y = 3*exp( -50*(x-.1)*(x-.1)) + 6*exp( -100*(x-.7)*(x-.7));break;
		case 5: y = 10*(pow(fabs(x-.2), 3) - pow(fabs(x-.3), 3) + 4*pow(fabs(x-.5), 3) - pow(fabs(x-.8), 3)); break;
		case 6: y = (fabs(x-.2) + fabs(x-.3) + fabs(x-.5) + fabs(x-.8)); break;
		case 7: y = sin((1./(x*x+.1)) + 1./(x-1.2)); break;
		case 8: y = 8*exp( -100*(x-.1)*(x-.1)) + 
					6*exp( -200*(x-.6)*(x-.6))+ 
					10*exp( -150*(x-.9)*(x-.9));break;
		case 9: y = 2*sin(8*exp( -100*(x-.1)*(x-.1))) + 
					6*exp( -200*(x-.6)*(x-.6))+ 
					10*exp( -150*(x-.9)*(x-.9));break;
	
	}
	return y;
}
